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Abstract. The mass spectrum of heavy pseudoscalar mesons, described as quark-antiquark 
bound systems, is considered within the Bethe-Salpeter formalism with momentum- 
dependent masses of the constituents. This dependence is found by solving the Schwinger- 
Dyson equation for quark propagators in rainbow-ladder approximation. Such an approx- 
imation is known to provide both a fast convergence of numerical methods and accurate 
results for lightest mesons. However, as the meson mass increases, the method becomes less 
stable and special attention must be devoted to details of numerical means of solving the 
corresponding equations. We focus on the pseudoscalar sector and show that our numerical 
scheme describes fairly accurately the ir, K, D, D s and r\ c ground states. Excited states are 
considered as well. Our calculations are directly related to the future physics programme at 
FAIR. 



1 Introduction 



The description of mesons as quark-antiquark bound states within the framework of the Bethe- 
Salpeter equation with momentum dependent quark masses, determined by the Schwinger-Dyson 
equation, is able to explain successfully many spectroscopic data [IHS]. Contrarily to traditional 
phenomenological models, like quark bag models, the presented formalism maintains impor- 
tant features of QCD, such as dynamical chiral symmetry breaking, dynamical quark dressing, 
requirements of the renormalization group theory etc. The main ingredients here are the full 
quark-gluon vertex function and the dressed gluon propagator, the calculation of which is en- 
tirely determined by the running coupling and the bare quark masses. In principle, if one were 
able to solve the Schwinger-Dyson equation within all pQCD orders, the approach would not 
depend on any free parameters. However, due to known technical problems, one restricts oneself 
to calculations of the few first terms of the perturbative series, usually up to the one-loop ap- 
proximation. The obtained results, which formally obey all the fundamental requirements of the 
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theory, are then considered as exact ones with, however, effective parameters. This is known as 
the rainbow-ladder approximation for the Schwinger-Dyson equation. The merit of the approach 
is that, once the effective parameters are fixed, the whole spectrum of known mesons is supposed 
to be described on the same footing: including also excited states. 

It should be noted that there exists other approaches based on the same physical ideas but 
not so sophisticated, e.g. employing simpler interactions, such as a separable interaction for 
the effective coupling [6|. Such approaches describe fairly well the properties of light mesons, 
nevertheless, investigation of heavier mesons and excited states, consisting even of light (u,d,s) 
quarks, requires implementations of more accurate numerical methods to solve the corresponding 
equations. Among other successful efforts in this realm the Refs. |10pil| must be also mentioned. 

In the present note we are going to apply the combined Schwinger-Dyson and Bethe-Salpeter 
(BS) formalisms to describe the meson mass spectrum including heavy mesons and excited 
states as well. Particular attention is paid to the charm sector which, together with the baryon 
spectroscopy, is a major topic in the FAIR research programme. Two large collaborations at 
FAIR [12^13] plan precision measurements. Note, that it becomes now possible to experimen- 
tally investigate not only the mass spectrum of the mentioned mesons, but also different processes 
of their decay, which are directly connected with fundamental QCD problems (e.g., £7(1) ax- 
ial anomaly, transition form factors etc.) and with the known problem of changing the meson 
characteristics at finite temperatures. The latter is crucial in understanding the di-lepton yields 
in nucleus-nucleus collisions at, e.g. HADES. All these circumstances require an adequate the- 
oretical foundation to describe the meson spectrum and the meson covariant wave functions 
(i.e. the BS partial amplitudes) needed in calculations of the relevant Feynman diagrams and 
observables. 

Our paper is organized as follows. In section [21 the parametrization of the gluon propagator 
is presented and the Schwinger-Dyson equation for the quark propagator is discussed. Section [3] 
deals with the Bethe-Salpeter equation. Numerical results are discussed in Section^ Conclusions 
are drawn in section O 



2 Propagators and Schwinger-Dyson equation 

The Bethe-Salpeter and Schwinger-Dyson equations in Minkowski space contain poles and 
branch-point singularities which strongly hinder the procedure of finding numerical solutions. 
Usually, to avoid these difficulties, one performs the Wick rotation and formulates the corre- 
sponding equations in Euclidean space, where all singularities in amplitudes and propagators 
are removed, so that the equations can be solved numerically. The known Mandelstam technique 
allows then to calculate matrix elements of observables which, being analytical functions of the 
relative energy, are the same in both Minkowski and Euclidean spaces. 

In our case we consider the Schwinger-Dyson equation for the quark propagator within 
pQCD with summing all diagrams up to one-loop. In calculations of diagrams the chiral sym- 
metry breaking is implemented from the very beginning [9]. The exact results, even only up 
to one- loop diagrams, after proper regularization and renormalization procedures are rather 
cumbersome for further numerical calculations. Nevertheless, in practical calculations one can 
employ reasonable parametrizations for the corresponding vertices and propagators to find solu- 
tions numerically. So, a phenomenological expression, inspired by calculations of the mentioned 
diagrams and preserving the requirements of the theory, for the combined running coupling and 
gluon propagator has been suggested in Ref. [9] 
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where the first term originates from the infrared (IR) part of the interaction determined by non- 
pertubative effects, while the second one ensures the correct ultraviolet asymptotics. Accordingly 
to the known fact that the contribution of the IR part is predominant for formation of bound 
states, in what follows we neglect the second term and restrict ourselves to the IR one. Then, 
as seen from Eq. ([1]), we are left with only two effective parameters, D and oj. 

The calculation of the renormalized Feynman diagrams leads to a fermion propagator de- 
pending on two additional parameters. In canonical calculations these are the renormalization 
constant Z 2 and the self-energy corrections £(p)- Usually, for further simplifications of calcu- 
lations, instead of Z 2 and 5J{p) one introduces other two quantities A{p) and B(p) in terms of 
which the quark propagator S q (p) reads as 

S-HP) = H ■ PA(P) + B( P ); S q (p) = '^Jffllfff ■ (2) 

Then with such a representation of the quark propagator the Schwinger-Dyson equation in 
Euclidean space has the form (cf. Refs. 019]) 

S q l {p) =il -p + m + ^ j ^4 [g 2 D(p-l)]^- lil S q {l)-i lJ , (3) 

where m is the bare quark mass and the effective kernel D(p — I) is 

D (*V = ^-'>'(^-^). (4) 

Note that ([3]) is a four dimensional integral equation. To solve it one usually decomposes the 
kernel over a complete set of basis functions, performs analyticaly some angular integrations and 
considers a new system of equations relative to such a partial decomposition. In our calculations 
we expand the interaction kernel into hyperspherical harmonics 

D((p - I) 2 ) = D(p 2 , 1 2 , \p\\l\ cos 7) = E D "(P> l )Cn(™s 7) , 

n 

1 

D n (p, I) = lj D(p 2 , l 2 ,p-l- tjC^tWl-^dt , (5) 



-1 



where C\{t) are the Chebyshev polynomials. For the employed ansatz of the gluon propagator 
the angular integration can be performed analytically leading to 



D n (p,l) = (n + l)8n 2 De- 



- — -I n+ i(z) - I n+2 (z) 



(6) 



D' n (p,l) = (n + l)An 2 ^ 1 -^- (7) 

with x = (p 2 + 1 2 )/uj 2 , z = 2pl/u 2 and I n (z) being the modified Bessel functions of the first 
kind. In Eq. ([7]) D' n (p,l) denotes the corresponding part of the kernel D(k 2 )/k 2 . Eventually, the 
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Figure 1. The solution of the Schwinger-Dyson equation for the quark mass m(p) = B(p)/A(p) as a function of 
momentum for different flavors, i.e. for the input bare quark masses m. At asymptotically large values of p the 
mass of the dressed quark tends to its bare value. 



system of equations to be solved reads as 



i 4 f dU A ( l ) r, < 

~ * f dllS B ( l ) n / A 

B{P) = m + *J^ PA*{l) + B*{l) DoM - 



The resulting system of equations on (JH]) is a system of one-dimensional integrals and can 
be solved numerically, e.g. by an iteration method. We found that iteration procedure for (J8j) 
converges rather fast and practically does not depend up on the choice of the trial functions for 
A{p) and B{p). The numerical solution of (J5|) with the effective parameters from Ref. is 
shown in Fig. [1] as the momentum dependence of the mass m(p) = B{p) / A{p) of the dressed 
quark for different quark flavors, i.e. for different bare masses: m = 0.005 GeV for u(d), fh = 0.115 
GeV for s and fh = 1 GeV for c quarks [7]. In principle, the asymptotic behavior of A(p) and 
B{p) for large momenta can be obtained directly from ([5]) and, as an additional check of the 
method, compared with the corresponding numerical results. 
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3 Bethe-Salpeter vertex function 

To determine the bound state mass of a quark-antiquark pair one needs to solve the Bethe- 
Salpeter equation. In rainbow-ladder approximation and in Euclidean space it reads [7l[9] 

r(p, P ) = (-i) J ^L^ S (k + )r(P, k)s(k^(g 2 D(p - k))„ v , (9) 

with P being the Bethe-Salpeter vertex function. The color structure has been factorized explic- 
itly. P means the total momentum of the bound state, p(k) is the relative momentum within 
the quark pair and k + = k + £P, = k + (£ — 1)P. The result does not depend on £ due to the 
covariance of the Bethe-Salpeter equation. We choose £ = 0.5. However, approximations might 
destroy covariance and one has to check a posteriori the stability of the results. 

Equation (jSJ) is written in matrix form, i.e. the vertex function P(P, k) is a 4x4 matrix 
and, therefore, may contain 16 different functions. The general structure of the vertex functions 
describing bound states of spinor particles has been investigated in detail, for example, in [14]. P 
can be expanded into functions which in turn are determined by angular momentum and parity 
of the corresponding meson: 

r(p*,p) = ^9a(p4,\P\)Ta(p). (10) 

a 

For pseudoscalar mesons there are 4 independent angular momentum functions T a {p) given by 

Tf \ 75 rr f x 7o75 rj* / \ 1 {Pi) m s 1 (pt) f ,,^ 

Vlo7r Vlo7r vlo7r |p| vlo7r \p\ 

Our choice differs from the standard expansion [9] . The main advantage of (jlip is orthogonality 
/ df2 p Tr[T m (p)T+(p)} = 5 mn . 

This property immediately allows to get a system of linear integral equations for the functions 
9a{PA,p)- As in the case of the Schwinger-Dyson equation, for these functions we employ also 
the hyperspherical harmonics 



which separates the angular dependence 

ft(p4,|p|)=X;fl/(p)X i)0(1) (xp) (13) 
j 

and, hence, allows to perform the angular integration leaving us with a set of coupled one- 
dimensional linear integral equations for the functions gf(p) 

f°° Hkk 3 . 

fl?(p) = / ^E4 m (f.%rw ( i4 ) 

" o 

where is the length of the Euclidean 4-momentum, p = \Jp\ + p 2 ■ The matrices A^p^p, k) 
can be calculated analytically. 

The series (|13p converges rather fast and in practice only few terms need to be taken into 
account. Then by choosing a suitable method of integration, e.g. Gaussian quadrature, the system 
of equations (|14p can be written as a system of homogeneous linear equations. Schematically, it 
can be written in the form 

9 = Kg, (15) 
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Figure 2. Integration contour of (|17[) . Branches of parabola are determined by Refc 2 - — —M 2 /4 — k 2 and 
Imfc± = ±Mfc. Since the integration domain is infinite, in concrete calculations one shall restrict to some large, 
but finite values of fc 2 ; the vertical line illustrates such a cut of the infinite domain in practical calculations. 



The condition det(K — 1) = is sufficient for the existence of a bound state. Hence, formally the 
zeros of the determinant det(-fT — 1) determine the solution of the BS equation, including also 
excited states. In our case, for the sake of increasing the accuracy, for the standard Gaussian 
quadrature we use a mapping 

k = k \±^- , (16) 
1 — x 

where ko weights the integrand and x £ [— 1, 1] is the integration variable (cf. Ref. [15J). 

Here an important moment is worth to be emphasized. Usually, when solving the BS equation 
for constituent particles, i.e. for particles with constant masses, the resulting system of partial 
equations is real. In case of momentum dependent masses the Bethe-Salpeter equation becomes 
complex and requires the knowledge of the quark propagator for complex momenta k± (where 
k± = \P ±k are the momenta of quarks) and, hence, requires to solve the Schwinger-Dyson 
equation also for complex momenta. For small meson masses (e.g. M ~ 500 MeV) only a small 
energy region contributes to the integral in Eq. (|14p . k < 1 GeV. In this case the propagator 
functions can be obtained by using the solution for real momenta, and afterwards A{p) and 
B(p) are calculated at complex momenta from ([8]) [3,7J. For heavier states, M > 1 GeV, the 
imaginary part of the quark momenta, Im k±, becomes rather large and the integrand in (JSj) 
rapidly oscillates as a function of k, hindering an accurate computation of the integral. 

To avoid this problem a witty trick has been suggested in Ref. [S]. It is based on the obser- 
vation that, as seen from the Schwinger-Dyson equation, all the values of the momenta k\ are 
located within a domain limited by a parabola, as shown in Fig. [2j 

Then, due to the Cauchy's theorem, it suffices to know the solution along this parabola to 
be able to compute it everywhere inside the corresponding domain^: 



27riA(zi) = j> dz 



Z — Z\ 



(17) 



Shifting the integration variable I — > p — I, ([8]) can be solved directly along this contour. For 
large real part of the momentum the asymptotic form of A and B can be also used. It should be 
noted that the described procedure of solving numerically the equations works quite well if the 
integration domain is chosen in a reasonable way, i.e. large enough to assure good asymptotics, 
but not too large, as the accuracy of the solution decreases with increasing integration domain 
(certainly, at a given Gaussian mesh). Note also, that asymptotic behaviors of the functions A 
and B is essential for heavier quarks, see Fig. [TJ 



x Many other methods may be employed to improve the accuracy of numerical calculations, see e.g. Ref. |16| 
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Figure 3. The bound state masses of a system qq x as a function of the bare mass rh x , the other quark being 
either a it quark with rh u — 0.005 GeV (the curve labeled as "itic") or a strange quark with fh s = 0.115 GeV (the 
curve labeled as "sx") or a charm quark with m c = 1.0 GeV (the curve labeled as "cx"). The effective parameters 
are ui — 0.5 GeV and D — 16 GeV -2 . The masses of pseudoscalar it K, D, D s and rj c mesons according to [17] 
are indicated at the right side. The vertical dashed lines mark the selected bare quark masses for it, s, and c. 



4 Numerical results 

As an example of our numerical study we exhibit in Fig. [3] the energy of the lowest bound states 
of a hypothetical meson qq x consisting of one given quark q with the mass known from the 
Schwinger-Dyson equation, bound with a second quark q x for which the input bare mass rh x 
is let to vary arbitrarily. The corresponding effective parameters have been chosen as uj = 0.5 
GeV and D = 16 GeV -2 [TIE] an d the bare masses for q correspond to (u, d, s) quarks, q=u 
(with fh u = 0.005 GeV), q=s (with rh s = 0.115 GeV) and q=c (with rh c = 1.0 GeV). This figure 
illustrates the whole mass spectrum of pseudoscalar mesons with masses up to 3 GeV/c 2 . So, 
if the q x quark corresponds to a c quark, then at the intersection of the vertical line rh x = m c 
(~ 1 GeV/c 2 ) with "ux" curve one obtains the D -meson (with the quark contents auc), with the 
"sx" curve the D s meson and with the " cx" curve - the rj c meson, respectively. It is worth noting 
that the "ux" curve crosses the rh u = 0.115 GeV line roughly at the same value of M qqx as the 
"sx" curve crosses the rh u = 0.005 GeV line, thus proving a check of consistency of the approach 
and, at the same time, describing correctly the lowest pseudoscalar us state corresponding to 
the K meson. It can be seen that even without a fine tuning of m SjC the meson mass spectrum 
is reproduced fairly well: 135 MeV (ir° meson), 497 MeV (K meson), 1870 MeV (D^ meson), 
1970 MeV (Df meson) and 2980 MeV (rj c meson). 

In nature, a pseudoscalar ss meson with simple quark structure does not exist. However, our 
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"sx" curve intersects with the bare mass of the s-quark, i.e. the approach predicts an existence of 
ss meson around 600 MeV. This can be considered as a "ghost" state, or as check of consistency 
of the approach. Apart from this circumstance, it is amazing that the simple two-parameter 
ansatz of the gluon propagator in the IR region together with the Schwinger-Dyson equation in 
rainbow approximation delivers such a quark propagator which, via the Bethe-Salpeter equation 
in rainbow-ladder approximation, results in such a nice description of pseudoscalar it, K, D, D s 
and T] c states. 

Clearly, a special parameterization of the gluon propagators with two adjusted quantities 
and the self consistent determination of three bare quark masses serve as input for obtaining 
finally five pseudoscalar meson masses. The straight forward extension to the scalar, vector 
and axialvector mesons without further adjustments will provide the ultimate test for the sub- 
tleties of the numerical implementation of the employed theoretical scheme and the inherent 
approximations. Work along this line is in progress. 

Masses of excited states, which correspond to radial excitations of the considered mesons, 
can be evaluated in a similar way, i.e. by finding of the next zeros of the determinant (|15p . As our 
analysis shows, the determinant changes monotonously with increasing mass. So that, naively, 
one would not expect a bound state of a uu system with mass larger than 800 MeV, since even at 
zero momenta, the maximum constituent quark mass is around 400 MeV, cf. Fig. (TJ However, in 
the deeply non-Euclidean domain, which corresponds to k\ ~ — M 2 /4 (Fig. [2]), the dynamical 
quark masses, contrubiting to the BS equation, increase till 600-700 MeV depending on M. 
Hence, the Schwinger-Dyson equation allows to understand the formation mechanism of such 
bound states which, from the constituent quark model point of view, can not be even predicted 
a priori. The mass of the first excited uu state is found to be 1080 MeV, i.e. significantly above 
the maximum mass from the Schwinger-Dyson equation alone. Analogously, for the cu system, 
the first excited state is found to be around 2530 MeV, which is in a good agreement with data. 
Similar results have been obtained in other groups, see Ref. [10J. 

Note that in solving the BS equation for the mass spectrum of mesons we obtain also the 
partial BS amplitudes which can be used in calculations of various dynamical observables, such 
as the meson life-time, transition form factors, the dependence of the meson widths up on 
temperature of an ambient medium etc. Such investigations are in progress and results will be 
reported elsewhere. 



5 Conclusion 



The method of solving the Schwinger-Dyson equation in rainbow-ladder approximation in Eu- 
clidean space by using the hyperspherical harmonics basis is proposed. The obtained numerical 
solutions are then used to solve the Bethe-Salpeter equation for the meson mass spectrum in 
a large interval of meson masses. In solving the Bethe-Salpeter equation a new set of basis 
functions has been used which allows a further easy decomposition of the Bethe-Salpeter vertex 
functions into hyperspherical harmonics basis 

The obtained mass spectrum for pseudoscalar mesons in a wide range, ranging from pions 
to r] c mesons, is in a good agreement with experimental data. Excited states were considered as 
well and found to be also in a good agreement with experimental data and with calculations by 
other groups. By solving the Bethe-Salpeter equation, the corresponding partial wave functions 
are also obtained which will allow, in future, to calculate a variety of observables related to 
physical programmes at, e.g. FAIR. 
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